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Abstract 

We consider numerical solution of Coupled Nonlinear Schrodinger Equation. We prove stabil- 
ity and convergence in the L2 space for an explicit scheme which estimations is used for implicit 
scheme and compare both method. As a test we compare numerical solutions of Manakov system 
with known analytical solitonic solutions and as example of general system - evolution of two 
impulses with different group velocity (model of pulses interaction in optic fibers). Last exam- 
ple, a rectangular pulse evolution, shows asymptotic behavior typical for Nonlinear Schrodinger 
Equation asymptotics with the same initial condition. 

1 Introduction 

A growing interest to encoded electromagnetic pulses propagation for long distances and processing 
is noted [TJ [2] ■ A lot of realistic models are based on the Coupled Nonlinear Schrodinger Equations 
(CNLSE): they are developed and used in many nonlinear optics problems such as polarization modes 
interaction [2l [3l |4] . For such modes some numerical schemes and codes are elaborated beginning 
from a celebrated integrable Ablowitz-Ladik one [51 . 

In papers [3 [S] IHl [13 [HI [12 authors investigate numerical methods for solving CNLS equation 
based on finite difference schemes. Such difference schemes are applied in ^13j to model vector spatial 
soliton behavior in a nonlinear waveguide arrays. 

There are numerical schemes which are unconditionally unstable, conditionally stable and un- 
conditionally stable. Generally could be said that we have two useful schemes: conditionally stable 
(Euler type) schemes and unconditionally stable (Crank-Nicolson type) schemes; we focus on these 
types of schemes. Most important thing is to use conservation laws for CNLS equations while the 
scheme constructed. There are a lot of possibilities to define conservation laws ^ but we choose in 
this paper the standard one given by [lOl [9] . 

Authors of [21 [Til [Zj compare explicit method with implicit one underlying that both of the 
methods have good and bad sides and there is no an uniform point of view. Conditionally stable 
methods do not need matrix formulation, but they need smaller space and time steps to assure 
stability. Unconditionally stable method need to solve system of equations, but bigger time step 
could be used. We would like to compare numerical solutions for implicit and explicit method under 
the scope of stability proof for explicit method. In paper [9] authors prove convergence and stability 
for NLS equation but do not analyze CNLS equations and how stability and convergence depends 
on initial amplitude (energy) of impulses. 

We consider a finite-difference scheme for the CNLSE which extends in a sense the results con- 
cerned linear Schrodinger equation [3], following the ideas successfully applied for the coupled 
Korteveg -de Vries (KdV) systems in [151 [H] ■ 
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Most important aim which we would like to achieve is to estimate stability and convergence 
parameters as a function of initial amplitude of pulses (energy) for the CNLS equations on a base of 
the norm originated from L2 space. Such result allows us to estimate time and space steps for the 
convergence regime inside which we could consider a solution as stable. 

The system of the CNLSE is considered in the form: 

idtU + iad^U + kd^^U + [a\U\^ + b\V\^] U = Q, (la) 

idtV - iad^V + kd^^V + [c\V\^ + d\U\'^] V = 0, (lb) 

where functions U and V are differentiable up to the second order and the following notations are 
used 

U 
V 

dtV = Ut 

^xx^ — Uxx 

and (for the example of two modes excited in the waveguide [I71[T5]) the parameter a is defined as 

iKi - Ki) • (3) 

This parameter describes the difference between group velocities of the modes. For simplicity of 
notation, but without losses of generality, we take a, 6, c,d,k,a > (in other case we should put 
everywhere in estimation of norm absolute value of this coefficients which could make equations 
more illegible). 

2 Conservation laws 

Let us prove first that the solutions of the system ([T]) satisfy the following conservation laws 

00 

J \U\'^dT = co72st, (4a) 

—00 
00 

J \V\^dT ^ const. (4b) 

— CSO 

If one writes the conjugate to the first of equations (H]) (equation for U amplitude) 

- id^U + iSdrU + kdrrU + [a\U\^ + b\V\'^]U = 0, (5) 
next multiply the result (O by —U and the equation ([T]) by U, it yields 

iUd^U - iSUdrU + kUdrrU + [a|C/p + b\V\^] UU = 0, (6a) 
iUd^U - iSUdrU - kUdrrU - [a\U\^ + b\V\'^] UU = 0, (6b) 



dU{x,t) 

= dt ' 

d'^U{x,t) 
dx-^ 



(2a) 
(2b) 

(2c) 
(2d) 
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with obvious relations 

UU = UU =\Uf. (7) 
As a next step, adding the equation ((6a)l to (j6bp 

ia^[/I7 - i5c»^;7[7 + /cc*^ [UdrU - UdrlJ] = 0, (8) 

and integrating the resuh, one arrives at: 

oo oo oo 

id^ J {UU)dT-i5 j dr{UU)dT + k j dr \UdrU - UdrU] dr = (9) 



+ CXD 



id^ / UUdT-iS\U\'^ + k [UdrU - UdrU] =0, (10) 



— OO 



Suppose the boundary conditions at both infinities 



lim t/ - 0, (11a) 

r— *±oo 

(lib) 



are imposed, 

oo 

id^ J |L/pdr = 0, (12) 



— OO 

we obtain the first conservation law in the form 

oo 



J \U\^dT = const. (13) 

— OO 

For V we built the conservation law in the same way 

oo 

J \V\^dT = const. (14) 

— oo 

Such conservation laws give us a possibility to define the basic norm in L2 space as 

00 

\\wf^ J {\U\^ + \V\^)dT^ const. (15) 
—00 

We would denote other norms by the same notation if it will not lead to a confusion. 

3 Discretization of the CNLSE system 
3.1 Explicit scheme 

Choose the explicit scheme in a simple form [19 with a second order discretization with respect to 
time and a third order one to space: 

- Ti^ iP - iP iP - 21 P + IP / \ 
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Let us derive conservation law for discrete CNLSE system (where U is complex conjugate of U). 

Multiplying equation (fTB]) by {Ul^^ + Ul) and complex conjugating of this equation by (C//^^ + C//). 
If we apply zero boundary conditions, we have 

N 



i=l 
N 



1=1 



N 

1=1 

N 
i=l 



3.2 Implicit scheme 

This six point scheme [10] 



n+l 



rn+l 



n+l 



rn+l 



rn+l 



4/l 



, Tin \Tjn 



0, 



(17a) 
(17b) 



(18) 
(19) 



is exactly Crank-Nicolson one [in]- This implicit scheme bases on the same finite difference as 
explicit scheme, we expect that implicit scheme converges to exactly solution quicker than explicit 
scheme. 

Elements in a node 1/2 are calculated by iterations. 



4 Stability 

We can separate real and imaginary part of U and V by put 

C/ = C + 
V = a + if3. 

Substituting this into ([T]) yields in four equation with real amplitudes 

+ (J^t + kT^^x + [a(e' + + + /?')] ?? = 0, 
-?7t - ar^t + Kxx + [a(C' + v") + + C = 0, 
at - aat + + [c(a2 + /J^) + ^(^2 + ry^)] /? = 0, 
-A + cr/3t + ka^^ + [c(a2 + /js-j ^ ^^^2 + ^2^j ^ ^ 

If we apply explicit scheme (jl6p to this equations, we can build time evolution matrix as 



(20a) 
(20b) 



(21a) 
(21b) 
(21c) 
(21d) 








rpj+l 

-^22 









\ 












rpj+l 

-'34 




rpj+l 

^43 


rpj+l 

J 44 


/ 



(22) 
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where 



Til — T22 — Si^r — 2^ {Si+l,r — <^i-l,r) , (23a) 



that acts in the vector space i?^, of the columns: 



V{c[(a^')' + (/?^')']+rf[(e 



r)2 



( ^{ \ 



V J 



(23c) 

(24) 
(25) 



Now we win proof a stabihty with respect to small perturbations of initial conditions dW [151 E] 

dW^+i = [|T'^'dW°. (26) 

k 

For the stability conditions we require that the operator Y[k must be bounded by a constant 
in a sense of Frobenius matrix norm 



< c. 



(27) 

(28) 
(29) 



For stability, sufficient condition could be write in form 

IIT'^'II <exp(pr), 

where p is a constant, but case with \p(t, h)\ < const < oo 

IIT'^II <exp(p(T,/i)T) 

is also sufficient for stability under condition t, /i — > and some dependence of t and h. 

In matrix T all elements are matrix with index of spatial grid. Now we could upper estimate p. 
first we upper estimate all of matrix elements using spectral norm of matrix [16j : 

Ark 

—r- + Ta m, 
h'' 

Ark 

— — + Ta m. 



IITllll 


< 


11^1211 


< 


117^21 II 


< 


117^2211 


< 


11^3311 


< 


11^3411 


< 


11^4311 


< 


117^4411 


< 



, („]\2l 



T(J 

Tcr 



4Tfc 

Ark 



)^1 



~h' 



(30a) 
(30b) 
(30c) 
(30d) 
(30e) 

(30f) 

(30g) 
(30h) 
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Taking into account and ([5D|) we can write 
TJ+^ < 4 + 4— 



2Tamax[(e^)2 + {^[f] + 2r6niax[(aJ")2 + {pr)^] 

i i 

+ 2rcmax[(a^)2 + + 2rdmax[(e)2 + (ry^X']. (3I) 

Taking into account the finite-difference conservation laws (jl7ap we write 

(32a) 
(32b) 



1=1 


i=l 


TV 
i=l 


AT 



now we can upper estimate norm (j3ip as 



|Tj+i 1 1 < 4 + 4— + — ^ + r2(a/„ + 6/„ + diu + ch))] < exp (pr). 



Finally we can write p as 



(7 16fc 

p = 4— + —5- + 4max(a, d)Iu + 4max(6, c)/^. 



(33) 



(34) 



This means that the fraction r/h must tend to (stability of this system of equation depends on 
initial energy quantity). For a better stability this condition requires that r — > much faster than 



5 Convergence 

In this section we proof that the terms at l.h.s of equations in finite difference form (jl6p converge to 
the equation ([T]), which solutions arc differentiable with respect to time and (twice) to x. We put 
into equation (jl6p a solution in the form 



= ES^ 



ri - 



dril 



eal + da\ 
\ e/3f + d(3i J 



(35) 



where ES is an exact solution and V is a difference between a numerical solution and the exact 
solution of CNLS equations. 

Below we show only one of the matrix component 



2h 



2h hT^ 



+ {a [(e^^' + d^if + [cTj^ + dry^y] + 6 [{eal + dc^ f + {ePj + d/?^'] } {ctjI + dr^j) = 0. (36) 
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Let us use the conservation laws p7ap for the equation ([T]) and the finite difference equation ([T 
that aUows us to write 



N 



N 



i=l 1=1 
N N 



(37a) 
(37b) 



or, in terms of the relation 



2h 



+ T {a [(e^D' + i^vir] + b \{ec4r + (e/Sf)'] } drjf 
- r {a [(eCD' + {erif) + b [{eajf + (e/J^)^] } drjf 



+ b 



+ T {a [(eef + rf^D" + (evi + dvff] + b [{eaj + daj)^ + (e/3f + dpff'^ } {eii + drji) 



(7- 



2h 



/l2 



+ r {a \{e^if + (er,^^] + b [(ea^)^ + (e/?^)^] } erjf 



(38) 



The right side of the equation for the differentiable solutions is of the order 0(r + ft, + /i^) [15], while 
to represent the left side we use the matrix T. Hence one arrives at 



dc 



i+1 _ tjh 



^V^ + r {a \{e^if + {e^ff] + b \{ealf + {e^if] } {evf + dri) 



■ {a \{e^ + d£,lf + {cT^l + d4f'\ + b ^eal + do^f + (epl + d/3|)2] } (e?/^' + dj^l) 

+ Q{t + h + h'^). (39) 



Let us define the norm of the numerical solution as 

\\U^\^h(^\Ul\^ . 
Now we upper estimate one element of the vector V-' 

Ue^^W < ||Ti+i||V^|| + T {alu + bh) II V II -T{alue + bhe) \W\\ + Q{t + h + h") 

< exp (iTp)||V"|| + r (a/„ + 6/„) HfZ-'H - r (a/„e + bhe)\\W\\ + Q{t + h + h^), 

Next we write the converge condition for all V-'+^ matrix components using the Frobenius norm 
up to the choice of the initial error \ \ = 0. 

\\V'+^\\ <Q + eiT + h + h^), (43) 

where i = 1, 2, 3, 4 and we define 

Q = 4|max(a,6,c,d)(/„ + /„)3/2 - max (a, 6, c, d)(/„e + /„e)^/^|- (44) 



(40) 



(41) 
(42) 
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6 Numerical results 



6.1 Nonlinear Schodinger equation 

For simple test we put c — d — and k = 0.5 in this case we have simple nonlinear Schrodinger 
equation. First we test NLS equation with initial condition U{0,x) — sech(a;) which pulse should 
not change shape during propagation (picture [T^) . 
As a full initial condition we set: 

U{0,x) = Aisech(a;), (45a) 
V{0,x) = 0, (45b) 

this lead to solve NLSE problem in form 

iUt + ^U^^ + \U\^U = 0. (46) 

When amplitude Ai — 1 the exact solution of equation pS)) is 

U{t, x) = Aisech{x) exp (it/2), (47) 




Figure 1: NLS equation with initial amplitude {Ai) a) 1 and b) 2 

Next we test NLS equation with amplitude A\ — 2\ the results are the same as in [2] (picture 
[T|d). For this case the exact solution of the equation (|46p takes more complicated form [5] 

[cosh(3a:) + exp(4jt) cosh(a;)] exp(ii/2) 
^'^ " cosh(4.T) + 4 cosh(2.T) + 3 cos(4t) ' ^ ' 

The numerical solution for the second case deviates more than the first one because it has bigger 
energy per pulse and we use for both cases the same time and space steps (see equation I34p. 



6.2 Manakov solitons 

Let us consider the soliton solutions of the Manakov system as examples to test the stability and the 
convergence of the explicit scheme. We start with the stability of the explicit scheme. At the picture 
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0,0018- 
0,0016- 




(a) (b) 



Figure 2: Error of numerical calculation for different time value a,) Ai = 1 and h) Ai — 2 

13] six cases of energy conservation behaviour are shown for the values of pr which are bigger than 
0.1: we observe very unstable results. When we decrease the time step, our solution is stabilized. 
It is very important to remark that p depends on initial amplitude of pulses. We choose the initial 
conditions as 

U{0,x) = Aisech{x), (49a) 
V{0,x) = A2sech{x). (49b) 

Table 1: Parameters for the numerical experiment "Manakov solitons" 



X 


-50 ... 50 


time 


15 


sigma 





a 


1 


b 


1 


c 


1 


d 


1 


Ai 


1 


A2 


1 


space steps 


1000 


time steps 


10000 - 1000000 



6.3 Collision of two solitons 

We make this experiment to compare the results with the " experiment 1" of the paper [TU] . General 
properties of such collisions of two solitons are well known 20J, hence we do not focus on details of 
this type of soliton interaction. 

6.4 Different group velocities 

The examples to be studied in this section are most important for us, because we interest in describing 
modes interaction in a nonlinear waveguide of different modes excited (different modes have different 
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41,0 
40,8 
40,6 
40,4 
40,2 
40,0 
39,8 
39,6 
39,4 
39,2 
39,0 





time step 


T*p 


A 


10000 


3 . 5 


B 


40000 


. 8 


C 


60000 


. 6 


D 


80000 


. 4 


E 


100000 


. 35 


F 


1000000 


. 035 



6 

time 




Figure 3: Stability of explicit method for Man- Figure 4: Manakov case:^i — A2 — I, nonlinear 
akov examples coefficients {a,b,c,d) = 1 

and (7 — 1. 





(a) \U\ 



(b) \V\ 



Figure 5: Inelastic collision of two solitons 
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group velocitys) [T71[TH] . This investigation could be useful not only for the situation described below, 
but for soliton trains interactions as well [21]. 

As initial conditions we took two sech-impulses and each equation have different nonlinear coeffi- 
cients. This could appear in a waveguide when two modes are excited with different group velocity. 
We set initial conditions for this case as: 

U{0, x) — Aisech(a; + Z^i) exp {ivix), (50a) 
V{0, x) = A2sech(a; + D2) exp {iv2x). (50b) 

Consider two solitons one with zero velocity and second with velocity greater than zero, we 
show this situations on picture [6l This two cases of pulses interaction differ only by group velocity 
(parameters for this pulses are in the table 16.4]) . 

This two impulses start from the same position (e.g. two modes excited in the waveguide). When 
velocity of second pulse is 0.7 (for given parameters) first impulses intercept part of energy thought 
nonlinear interactions and move with second pulse with average velocity for both pulses (picture 
|6K)- When the velocity of second impulses is high enough we have two pulses which moves with 
different group velocity see picture [GId (nonlinear interaction between pulses happen only at the start 
of propagation). 



Table 2: Parameters for numerical experiment "Different group velocity' 





(a) 


(b) 


X 


-30 ... 30 




time 


40 


40 


cr 








a 


1 


1 


b 


1/3 


1/3 


c 


1 


1 


d 


1/3 


1/3 




1.2 


1.2 


A2 


1.4 


1.4 


Vl 


0.7 


0.95 


V2 








h 


0.2 


0.2 


T 


0.02 


0.02 



6.5 Difference between explicit and implicit scheme 

If we would like to have more space details (smaller h), in explicit scheme we must take adequate 
smaller time step to assure stability of numerical scheme. This have very big influence on calculation 
time. In an implicit scheme we could use an iteration method [191 [7| and shorter time of calculation 
is achieved. On the picture [7| we show numerical results for both methods. Note, the time step for 
the explicit scheme is much larger than for the implicit scheme but each step in the implicit scheme 
needs 2-4 iterations. 

6.6 Rectangular pulse decay 

In this subsection we engaged in asymptotic solution of Nonlinear Schrodinger Equation. As initial 
condition we choose rectangular pulse as in [22] (see picture [Hd). 
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(a) (b) 
Figure 6: Two cases of impulses with different group velocity (for parameters see table 
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X 



(a) 



(b) 




'AAav 



-30 -20 -10 10 20 30 

X 



(c) 



(d) 



Figure 8: Rectangular pulse evolution. On picture a) there is a 3D plot on picture b),c) and d) there 
are intersection of plot a) in different time. 
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Table 3: Parameters for numerical experiment for implicit and explicit method 



Parameter 


explicit 


implicit 


space step 


300 


300 


time step 


1000000 


2000 


time 


40 


40 


X 


-30. ..30 


-30. ..30 


Ai 


1.5 


1.5 


A2 


1.5 


1.5 


a 


1 


1 


b 


0.2 


0.2 


c 


1 


1 


d 


1.6 


1.6 


a 


0.3 


0.3 



We test our method for this kind of initial pulse. First we made two calculations with different 
time steps (for the first t — 2e^'* and r = 2e~^ for the second one) and compare errors between 
(picture [9]i). The difference between calculations is of the order determined by rp. Next we check 
the energy conservation (see picture [QId). 



% 0,000 



-30 -20 



20 30 



(a) 



,0 6,OOE-011 



4,OOE-011 



fl) 2,OOE-011 




(b) 



Figure 9: a) Difference between two numerical solutions with different time steps (t = 2e ^ and 
T = 2e^^. b) The energy deviation (to the initial state) for the time step r = 2e"'''. 

For the presented time (t — 4) we expect that the transient stage from rectangular shape of the 
pulse to the asymptotic behaviour of solution is realized. On the picture [TO] the evolution (maximum 
amplitude) of rectangular pulses with different width is presented. If the pulse have smaller width 
then the decay is faster. 



7 Conclusion 

We study a convergence and stability conditions which would be useful to estimate other numerical 
schemes for NLS as well. Additionally these conditions could be compared with ones for numerical 
methods for other nonlinear equations [TO] and show a character of CNLS equations (from numerical 
side). Most important results which we show in this paper is how to estimate a time step for the 
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Figure 10: Maximum amplitude for different pulse width (in round brackets) for time step t = 2e 
(maximum time in dimensionless unit is equal 4). 



explicit method due to initial energy of pulses. The results could be adapted to other methods based 
on this basic scheme which we use [8] . There is a possibility to enhance this method to third or high 
order in time and this method could be used as a starting point (to calculate lacked points of grid). 
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